Disentangling the causal relationship between rabbit growth and cecal microbiota through structural equation models

Background The effect of the cecal microbiome on growth of rabbits that were fed under different regimes has been studied previously. However, the term “effect” carries a causal meaning that can be confounded because of potential genetic associations between the microbiome and production traits. Structural equation models (SEM) can help disentangle such a complex interplay by decomposing the effect on a production trait into direct host genetics effects and indirect host genetic effects that are exerted through microbiota effects. These indirect effects can be estimated via structural coefficients that measure the effect of the microbiota on growth while the effects of the host genetics are kept constant. In this study, we applied the SEM approach to infer causal relationships between the cecal microbiota and growth of rabbits fed under ad libitum (ADGAL) or restricted feeding (ADGR). Results We identified structural coefficients that are statistically different from 0 for 138 of the 946 operational taxonomic units (OTU) analyzed. However, only 15 and 38 of these 138 OTU had an effect greater than 0.2 phenotypic standard deviations (SD) on ADGAL and ADGR, respectively. Many of these OTU had a negative effect on both traits. The largest effects on ADGR were exerted by an OTU that is taxonomically assigned to the Desulfovibrio genus (− 1.929 g/d, CSS-normalized OTU units) and by an OTU that belongs to the Ruminococcaceae family (1.859 g/d, CSS-normalized OTU units). For ADGAL, the largest effect was from OTU that belong to the S24-7 family (− 1.907 g/d, CSS-normalized OTU units). In general, OTU that had a substantial effect had low to moderate estimates of heritability. Conclusions Disentangling how direct and indirect effects act on production traits is relevant to fully describe the processes of mediation but also to understand how these traits change before considering the application of an external intervention aimed at changing a given microbial composition by blocking/promoting the presence of a particular microorganism. Supplementary Information The online version contains supplementary material available at 10.1186/s12711-022-00770-2.

conducted by Sánchez et al. [7] identified different quantitative trait loci (QTL) for both these traits, supporting the idea that the genetic control of the growth of animals fed under a restricted regime differs from that of animals fed ad libitum. This result is also supported by the estimate of the genetic correlation of 0.59 between these traits that was reported by Piles and Sánchez [5]. In line with results of previous research in human [8] and livestock [9][10][11] populations, a recent study conducted by Velasco-Galilea et al. [12] in the same rabbit population as used in [7] suggested that, in rabbit, cecal microbial composition and diversity are heritable traits, and identified several host genomic regions that affect cecal microbial composition. However, since most of these studies are based on associations or correlations, causal relationships between the microbiome and the traits cannot be implied. Thus, which fractions of the genetic effect on a phenotypic trait (such as ADG AL and ADG R ) are exerted directly versus indirectly, e.g. by promoting a specific microbial profile in the gut that has a favourable or unfavourable effect on the trait, remains unknown. As Gardiner et al. [3] pointed out in their review: "some specific examples are outlined where an increased abundance of certain microbiota members may be a result of improved productivity traits rather than the cause of them". In addition, when analysing the impact of the microbiome on a host phenotype, this effect can be confounded with host genetics effects if the genetic association between the microbiome and the phenotype is not considered in the model. How important these indirect host genetics effects that are exerted through the gut microbiota are relative to the direct genetic effects on growth and FE has not been evaluated in rabbits. Direct and indirect effects could be opposite in sign, leading to a null effect of the genotype on the trait if they are of similar magnitude. This could be the case of a microorganism that has a negative effect on a phenotypic trait and for which its relative abundance is associated with a host genetic marker that is linked to a gene that positively affects the same trait through some metabolic process. If both effects are of similar magnitude, the effect of that marker would not be captured by a standard GWAS for the trait. Structural equation models (SEM) [13] can disentangle the complex interplay that exists between the host genome and its microbiome. SEM can also help to predict the effect that an external intervention would have on the system, such as one directed at inhibiting or blocking the effect of certain microorganisms. Multiple-trait animal models (MTAM) capture correlations or associations among traits but do not inform about their causal relationship. In contrast, SEM allow the decomposition of the total effect (genetic and environmental) on a trait into direct and indirect contributions. Thus, SEM can provide better insights in the relationships and the biological mechanisms among traits than MTAM.
The objectives of this research were: (1) to identify the operational taxonomic units (OTU) with abundances in the fecal microbiome that impact rabbit growth (ADG AL and ADG R ), (2) to quantify the relative importance of direct and microbial-mediated host genetic effects on ADG AL and ADG R , (3) to estimate the heritability of the abundance of the most important OTU, and (4) to assess the importance of direct genetic relative to total genetic covariances of ADG AL and ADG R with abundance of OTU.

Animals and phenotypic data
The 412 animals included in the present study belong to a line of rabbits that has been selected for post-weaning growth since 1983 and that is commonly used as a terminal sire line in the three-way crossbreeding scheme for rabbit meat production [14]. They were randomly selected from five batches of a larger experiment [5]. Most were produced in four batches in a semi-open-air facility of the Institute of Agrifood Research and Technology (IRTA, Barcelona Spain) during the first semester of 2014, and the remaining were produced in a single batch in another facility of IRTA under controlled environmental conditions in the spring of 2016. The animals bred in the first facility were housed in groups of eight kits from weaning (32 days of age) until the end of the fattening period (66 days of age). The kits raised in the second facility were housed in groups of six kits and their growing period was slightly shorter (from 32 to 60 days of age). A maximum of two kits from the same litter were assigned to the same cage to avoid confounding between cage and maternal effects. Apart from the afore-mentioned differences, all animals were raised under the same management conditions and were fed with a standard pelleted diet supplemented with antibiotics, except for 23 kits raised in 2016 that received the same food but free of antibiotics. Water was provided ad libitum and feed was supplied once per day in a feeder with three places during the 4-to 5-week fattening period. At weaning, the animals were randomly assigned to the ad libitum (AL) or restricted (R) feeding regime. The amount of feed supplied to the animals under R during each week for each batch was computed as 0.75 times the average FI of kits under AL from the same batch during the previous week, plus 10% to account for the increase in feed intake as the animals grow [15]. To generate homogeneous groups, kits under each feeding regime were categorized into two groups according to their individual weaning weight (WW) being greater or smaller than 700 g.
Individual body weight (BW) was recorded weekly for all animals. Individual ADG was computed as the slope of the within-animal regression of BW on age. This trait was computed for individuals under each feeding regime, thus obtaining ADG under AL (ADG AL ) and R (ADG R ). The number of animals and descriptive statistics for both traits are in Table 1.

Cecal sampling, microbial DNA extraction, and bioinformatic processing
For each animal, a cecal sample was collected in a sterile tube immediately after slaughter, kept cold in the laboratory (4 °C) and, stored at -80 °C until DNA extraction. Full details of the DNA extraction, amplification, library preparation and sequencing are given in previous reports [16,17]. Briefly, DNA integrity/purity was checked according to the protocol of Desjardins and Conklin [18], and then DNA was amplified following Parada et al. [19].
Raw reads were processed with the QIIME software version 1.9.0 (https:// github. com/ bioco re/ qiime/ relea ses/ tag/1. 9.0) [20], as described in [16]. Contigs with a quality score smaller than Q19 were removed. The UCHIME algorithm [21] was used to detect and remove chimeric sequences. Filtered sequences were clustered into OTU and those that were detected in less than 5% of the samples and with a count sum less than 0.01% were discarded. The remaining 946 OTU were normalized using the cumulative sum scaling (CSS) method [22]. Taxonomic assignment of each OTU was obtained by mapping the sequences against the Greengenes reference database. Raw sequence data were deposited in the sequence read archive of NCBI under BioProject accession number PRJNA524130.

Host DNA extraction and single nucleotide polymorhism (SNP) genotyping
Using the MN Nucleospin Tissue kit (Macherey-Nagel, Germany), host DNA was extracted from liver samples collected at slaughter. DNA integrity and purity were measured following the protocol of Desjardins and Conklin [18]. The Affymetrix Axiom OrcunSNP array (199,692 SNPs) was used to genotype 412 rabbits. Quality control of the SNPs was performed with the PLINK software (version 1.9) [23] according to four criteria: (i) individual call rate > 0.90; (ii) SNP call rate > 0.95; (iii) SNP minor allele frequency (MAF) > 0.05; (iv) and only autosomal SNPs with known positions in the OryCun2.0 assembly [24] were used. The final dataset contained genotypes for 114,604 SNPs on 412 rabbits.

Statistical analyses
With the main objective of estimating the effect of the microbiota on growth, while taking into account that these two traits are genetically correlated, the SEM represented in Fig. 1 was implemented to assess the causal relationship of each OTU (M) with ADG AL and ADG R . In this model, the host genotype directly affects the phenotypes (G → ADG AL and G → ADG R ) but also the gut microbiome (G → M), while, the microbiota also affects the phenotypes (G → M → ADG AL and G → M → ADG R ). In Fig. 1, AL←M and R←M are the structural coefficients and indicate the effect of M on ADG AL and ADG R , respectively, while keeping G constant. Environmental factors (E) affect M and the phenotypes.
For each OTU, the SEM was implemented using a Bayesian analysis of three phenotypes: ADG AL , ADG R , and OTU (referred as M). The analysis was repeated for each OTU, i.e. 946 times. The SEM for ADG AL , ADG R and the m th OTU (m = 1, 2, …, 946) for the n animals can be written as: where y is the vector of phenotypic records containing observations for ADG AL , ADG R , and M (with ADG R missing for animals on AL and ADG AL missing for (1) y = ( ⊗ I n )y + Xb + Z l l + Z c c + Z u u + e,  where G is the genomic relationship matrix and L 0 , C 0 , G 0 , and R 0 are symmetric covariance matrices for litter, cage, genomic and residual effects, respectively, with elements equal to (only the upper triangular matrix is shown):  . Note that, the parameters of the SEM are identifiable for any acyclic structure among traits, because R 0 is diagonal [25]. Solving Eq. (1), the SEM becomes: where I 3n is a (3n × 3n) identity matrix. Thus, the reduced model is: which corresponds to a MTAM with: and and In the SEM, u AL , u R, and u M (component vectors of u) represent the genomic effects that directly affect ADG AL , ADG R , and M, respectively. They can also be described as the effect of the genome on a trait, while holding the value for the remaining traits constant, i.e. the direct effects of the genome on a trait free from genetic effects that are mediated through other phenotypic traits that have a causal influence on the target trait [26]. In turn, AL←M u M and R←M u M represent the indirect genomic effects that are mediated by M on ADG AL and ADG R , respectively, such that u * AL = AL←M u M + u AL and u * R = R←M u M + u R are the overall genetic effects exerted by the genome of an individual on ADG AL and ADG R , respectively, through all the paths. Thus, the indirect effects on ADG AL and ADG R correspond to the direct effects of genes exerted over M, which, in turn, influence ADG AL and ADG R .
Components of G * 0 can be written as: where σ 2 u * ,AL , σ 2 u * ,R , and σ 2 u * ,M are the total genetic variance of ADG AL , ADG R , and M, respectively; σ 2 u,AL , σ 2 u,R , and σ 2 u,M are the genetic variance of direct effects on ADG AL , ADG R , and M, respectively; and σ u,AL,M and σ u,R,M are the covariances between direct effects, which represent the effects of genes that directly affect both traits (i.e. pleitropic ADG AL and ADG R effects) or of genes that affect one of the two traits and are in linkage disequilibrium with each other, respectively. Structural coefficients AL←M and R←M were as defined above. And and are the total covariance between ADG AL and ADG R , between ADG AL and M, and between ADG R and M, respectively. Note that if AL←M × σ 2 u,M and R←M × σ 2 u,M have opposite signs but equal magnitude to σ u,AL,M and σ u,R,M , respectively, the total genetic covariances between the ADG traits and M could be null. However, σ u,AL,R + AL←M × σ u,M,R + R←M × σ u,M,AL + R←M × AL←M ×σ 2 u,M is the indirect genetic association between ADG AL and ADG R (i.e., the part of the covariance that is due to microbial-mediated effects). The same Eqs. (14), (15), and (16) for the partition of the total genetic variances and covariances also apply to the other variance components of the model.
The SEM was implemented under a Bayesian approach via Gibbs sampling using the Gibbsf90 software of the blupf90 family of programs [27]. Single sampling processes of 1,000,000 iterations were run for all the models, discarding the first 300,000 iterations of each chain and saving 1 of every 100 samples. Posterior marginal inferences on variance components were performed with the Postgibbsf90 software [27]. Structural coefficients were considered to be statistically different from 0 when the high posterior density interval at 95% of confidence level (HPD95%) did not include the 0.

Recursive effects
Our SEM approach allowed the identification of 138 OTU with structural coefficients that were statistically different from 0: 82 OTU had an effect on ADG AL , 80 on ADG R , and 24 on both traits. Of these 138 OTU, 104 had a positive effect, with 45 affecting ADG AL , 45 affecting (14) ADG R , and 14 affecting both traits. Note that the resolution of the 16S rDNA locus that was used in this study to characterize the cecal microbiota allowed annotation of only about half of all sequences at the family level. Thus, most of the OTU with a positive effect on ADG could not be assigned at the family level (57, 42, and 50% for ADG AL , ADG R , and both traits, respectively). In the three cases, the majority of these OTU belong to the Lachnospiraceae and Ruminococcaceae families. For these two families, 22 and 17%, respectively, of the total number of OTU had a positive effect on ADG AL , 22 and 22% on ADG R , and 29 and 14% on both traits (Fig. 2).
The marginal posterior means of the positive structural coefficients for the effect of M (i.e. the specific OTU) on ADG AL ranged from 0.369 to 1.271 (g/d CSS-normalized OTU units) (see Additional file 1: Table S1), and from 0.452 to 1.859 (g/d CSS-normalized OTU units) for the effect of M on ADG R (see Additional file 1: Table S2). In turn, the marginal posterior means of the negative structural coefficients for the effect of M on ADG AL ranged from − 0.770 to − 1.907 (g/d CSS-normalized OTU units) (see Additional file 1: Table S1), and from − 0.707 to − 1.929 (g/d CSS-normalized OTU units) for the effect of M on ADG R (see Additional file 1: Table S2). To identify and compare OTU with a substantial effect on ADG AL and ADG R , the structural coefficient estimate for each OTU and trait was expressed as the expected change in number of phenotypic standard deviations (SD) for ADG for a one-unit increase in CSS-normalized OTU units by dividing the structural coefficient by the SD of the trait (Fig. 4).
Given the large number of OTU with structural coefficients that were statistically different from 0, only those with an absolute effect equal or higher than 0.2 SD units are described in detail. Estimates of these OTU, which we will refer to as relevant OTU, are listed in Table 2 (for ADG AL ) and 3 (for ADG R ). Among the OTU with a negative effect on the production traits, the New.ReferenceOTU4438, a member of the S24-7 family, was the OTU with the largest effect on ADG AL (− 1.907 CSS-normalized OTU units) ( Table 2) and the New.ReferenceOTU369, taxonomically assigned to the Desulfovibrio genus, was the OTU with the largest effect on ADG R (− 1.929, g/d, CSS-normalized OTU units) ( Table 3). Among the OTU with a positive effect on ADG R , a member of the Ruminococcaceaae family (New.ReferenceOTU1337) had the largest effect, with an estimated increment of 1.8 g/d per CSS-normalized OTU unit. Of the relevant OTU with a positive effect on ADG AL , OTU 641783 that had the largest effect and also belonged to the Ruminococcaceaae family, with an     for ADG AL and ADG R , respectively) and indirect effects ( 2 × AL←M × σ u,M,AL + 2 AL←M × σ 2 u,M and 2 × R←M × σ u,M,R + 2 R←M × σ 2 u,M for ADG AL and ADG R , respectively) for the most relevant OTU regarding their effects on ADG AL (Fig. 5) and ADG R (Fig. 6). In general, the distributions of the total and the direct genetic variances were nearly the same because of the small contribution of indirect effects. However, two relevant OTU of the Clostridiales order (New.Refer-enceOTU1337 and New.ReferenceOTU381 belonging Table 3 Posterior means, posterior medians, and 95% highest posterior density intervals (HPD 95% ) of structural coefficients (in g/d, CSS-normalized OTU units) for the relevant OTU on average daily gain of rabbits fed on restriction, along with their posterior means expressed in SD units and their assignment at the lowest taxonomic level to the Ruminococcaceae and Lachnospiraceae families, respectively) exerted a negative indirect effect on ADG R (Fig. 6), while OTU 209947, also of the Clostridiales order, negatively affected ADG AL . OTU that belong to other orders also affected ADG AL (New.Referen-ceOTU4438, Bacteroidales order) and ADG R (New. ReferenceOTU369, Desulfovibrionales order), but the contributions of their indirect effects to the total genetic variance were smaller.

Heritability estimates for ADG and the relevant OTU
Average heritability estimates across all 946 analyses for the two ADG traits were low to moderate, with posterior means (SD) of 0.215 (0.104) and 0.157 (0.101) for ADG AL and ADG R , respectively. Means (SD) of the marginal Fig. 5 Estimates of the total additive genetic variance ( σ 2 u * ,AL ) for average daily gain of rabbits fed ad libitum (ADG AL ), and of direct ( σ 2 u,AL ) and indirect genetic effects ( 2 × AL←M × σ u,M,AL + 2 AL←M × σ 2 u,M ) that contribute to this variance for the most relevant OTU Fig. 6 Estimates of the total additive genetic variance ( σ 2 u * ,R ) of the average daily gain of rabbits under restricted feeding (ADG R ), and of direct effects ( σ 2 u,R ) and indirect effects ( 2 × R←M × σ u,M,R + 2 R←M × σ 2 u,M ) that contribute to this variance of the most relevant OTU posterior distributions for the heritabilities for the most relevant OTU that affected ADG AL , and ADG R , or both, are shown in Fig. 7. In general, the estimates of these heritabilities were also low to moderate. Posterior mean heritabilities for the relevant OTU for ADG AL ranged from 0.083 for OTU 522353 to 0.222 for OTU 209947, both members of the Clostridiales order. Posterior mean heritability estimates for the relevant OTU for ADG R ranged from 0.046 for a member of the Phascolarctobacterium genus (New.ReferenceOTU570) to 0.234 for the New. ReferenceOTU369, which belongs to the Desulfovibrio genus. The most heritable families were Desulfovibrionaceae (0.23) and Methanobacteriaceae (0.24), while the least heritable families were Veillonellaceae (0.04) and Peptococcaceae (0.08).

Direct genetic covariances between ADG and OTU
Marginal posterior distributions of the total additive genetic covariance between ADG AL and M ( σ u * ,AL,M ), of the covariance due to direct effects ( σ u,AL,M ), and of the covariance due to indirect effects of M on ADG AL ( AL←M × σ 2 u,M ) are represented in Fig. 8. The same distributions are shown in Fig. 9 for ADG R ( σ u * ,R,M , σ u,R,M and R←M × σ 2 u,M , for total, direct, and indirect covariances, respectively). Total and direct additive genetic covariances were very similar for the New.Referen-ceOTU4568, OTU 339336, and New.ReferenceOTU1863 (all members of the Clostridiales order), due to the almost zero indirect effects. The indirect effects of the OTU 209947, which belongs to the Clostridiales order, on ADG AL were substantial and negative, while the indirect effects of members of the Desulfovibrio genus (New. ReferenceOTU369) and the Lachnospiraceae family (New.ReferenceOTU381) on ADG R were substantial and negative and positive, respectively.

Discussion
In this study, SEM were implemented to assess causal relationships between the cecal microbiota and two traits related to rabbit growth and FE. The effect that a specific OTU exerts on the studied phenotypes represents the indirect effects (i.e. the effect of the OTU on the trait while holding host genetics and environmental effects constant). It is measured by the structural coefficient, which quantifies the expected change in the phenotype by a one-unit increase in CSS-normalized OTU units of a specific OTU. When direct and indirect effects are of the same magnitude but have opposite signs, the total effect on the trait becomes null. In such a scenario, a GWAS can fail to identify the genomic regions that affect the production trait through multiple paths. This hypothesis was tested in a previous study conducted by Tiezzi et al. [28], who evaluated how host genetics can affect fat deposition in pigs by affecting gut microbiome composition, which then results in a change in the phenotype. In that study, it was demonstrated that the genes that do not affect the phenotype directly can be identified in GWAS in which effects mediated by the microbiome were included in the model. This highlights the importance of knowing the direct and indirect effects on a phenotype when designing sustainable strategies to improve productivity and sustainability in livestock. To our knowledge, the importance of direct versus indirect effects on growth rate through microbial diversity and composition have not yet been quantified. The only related study was by Saborio et al. [29], who implemented a SEM approach to estimate the effect of the relative abundance of rumen microbes on methane production in dairy cattle.
Most of the OTU that had statistically significant structural coefficients and that positively affected ADG AL and ADG R belong to the Lachnospiraceae or Ruminococcaceae families. Conversely, the OTU that had statistically significant structural coefficients and that negatively affected ADG AL belong to a wider set of families: Ruminococcaceae, Lachnospiraceae, Methanobacteriaceae, and Erysipelotrichaceae. Lachnospiraceae and Ruminococcaceae are both members of the Clostridiales order. Several members of this order have previously been found to be associated with FE across different livestock species, including pigs [30,31], chickens [32] and beef cattle [33]. Erysipelotrichaceae is a family that is associated with lipid metabolism and has been linked to inflammation [34]. However, its increase in abundance has been associated with the most but also with the least efficient pigs, depending on the sampling origin that each study chose for gut microbiota characterization [35].
As previously mentioned, given the large number of OTU with statistically significant structural coefficients that were identified by our SEM approach, our discussion will be focused on the OTU that had a substantial effect (equal to or larger than 0.2 units of SD of the trait) on the phenotype, i.e. the relevant OTU. We identified 15 and 38 OTU with relevant effects on ADG AL and ADG R , respectively. The estimated effects for the relevant OTU for ADG R were larger than those for the OTU affecting ADG AL . Thus, only 18% of the OTU that affected ADG AL had a structural coefficient larger than 0.2 SD units for ADG AL , while 47% of the OTU that affected ADG R met this criterion. This could indicate that growth rate is mainly determined by host genetics whereas FE is determined by both the host genetics and the gut microbiota. This is in agreement with results reported by Velasco-Galilea et al. [36], who found that the predictive value of cecal microbiota was higher for ADG R than for ADG AL .
It is worth highlighting the effects of some other OTU on ADG AL or ADG R . For instance, a member of the Phascolarctobacterium genus (New.ReferenceOTU570) had a substantial effect on ADG R (− 1.70 [− 3.194, − 0.475]). Some species of this genus produce acetate and propionate short-chain fatty acids (SCFA), which act as energy sources. Propionate is a gluconeogenic substrate in the liver and intestine, while acetate contributes to the synthesis of lipids [37,38] and has been shown to play a direct role in the regulation of appetite in mice [39]. Another OTU with an abundance that had a substantial positive effect on FE was the New.ReferenceOTU381 (Lachnospiraceae family). Velasco-Galilea et al. [36] reported five OTU of the Lachnospiraceae family (including the New.ReferenceOTU381) to be among the most informative OTU to predict rabbit growth and FE. Two and three of these OTU were, respectively, positively and negatively correlated with average daily residual feed intake in rabbits fed ad libitum.
Genetic variances for ADG AL and ADG R in the SEM can be interpreted as the variance due to direct genetic effects on the traits free from the genetic effects mediated by M that have a causal effect on them. Indirect genetic effects can act by reducing the total genetic variance if they are strong enough and opposite in sign to the genetic covariance between M and the phenotype. This could make it difficult to assess the genetic determinism of a trait affected by M based on MTAM. However, a SEM enables the different sources of genetic variation to be disentangled. In addition, it enables prediction of the effects of external interventions on a set of variables. In our study, an external intervention could involve promoting or blocking the presence of some cecal microbe.
One member of the Ruminococcaceae family (New. ReferenceOTU1337) and one member of the Lachnospiraceae family (New.Reference381) had relevant effects on the total genetic variance for ADG R , with indirect effects that reduced the total variance, although the estimates of the structural coefficients for these OTU were positive. This was due to the negative genetic correlation between ADG R and the abundance of these OTU. It is worth mentioning that some members of the Lachnospiraceae family are butyrate-producing bacteria that have a beneficial effect on animal gut health [40] and have been previously identified to have a strong effect on FE traits in pigs [41].
Estimates of total heritability for ADG AL and ADG R are similar to those estimated in previous studies for the same rabbit population, with a posterior mean of 0.21 and 0.08, respectively, by Piles and Sánchez [5] and with a posterior mean of 0.15 and 0.09, respectively, by Velasco-Galilea et al. [36]. Conversely, estimates of direct heritability based on the SEM were slightly greater than estimates of total heritability for both traits. This could be due to the negative effects of host genetics exerted through the microbiota on these growth traits. In general, estimated heritabilities for the OTU with the largest structural coefficients were low to moderate, suggesting a weak host genetic control of the rabbit cecal microbiota, which is in agreement with results of Velasco-Galilea et al. [12], who assessed the host genetics effects on the rabbit cecal microbiota by means of Bayes factors. The highest heritability [0.234 (0.091)] was estimated for New.ReferenceOTU369, which belongs to the Desulfovibrio genus and is one of the OTU that had a relevant effect on ADG R . This OTU was previously reported as heritable by Velasco-Galilea et al. [12].
In SEM, the genetic covariance reflects the association between direct effects and can be considered to be due to genes that directly affect two traits simultaneously or to linkage disequilibrium between pairs of genes that each affect one of the two traits [13]. However, there is a second source of association between the host phenotype and M because the host genetic effect on M also affects the host phenotype indirectly. This second source of covariation could even have an opposite sign compared to the first source (i.e., the covariance between direct genetic effects). A member of the Desulfovibrio genus (New.ReferenceOTU369) impaired ADG R , and the total indirect effects that this OTU was estimated to exert on the phenotype acted by reducing the total effects. In this particular case, direct and total covariances were opposite in sign. Thus, the genetic correlation between this particular OTU and ADG R would be positive if only the effects of host genetics were considered. On the contrary, when the total effects are considered, the relevant negative effect of the Desulfovibrio genus on ADG R can be captured by the negative genetic correlation between the OTU and growth. This result is in agreement with previous reports in pigs [30,42] in which the abundances of members of the Desulfovibrio genus were estimated to be negatively correlated with FE.

Conclusions
Our study highlights the importance of knowing the direct effects that host genetics exerts on a phenotype, as well as the indirect effects that host genetics exerts through the gut microbiota, not only to fully describe the processes of mediation, but also to understand the change in a phenotype that can result from a modification of the microbial composition through an external intervention (e.g., by blocking/promoting the presence of a particular microorganism). This is the first study to evaluate the direct and indirect effects exerted by host genetics on growth. Our results show that 63% of the OTU with abundances that had relevant effects on ADG R had positive effects on this trait. Abundance of one member of the Desulfovibrio genus exerted the largest negative effect on ADG R , followed by abundance of one member of the Ruminococcaceae family, which positively affected this trait. In contrast, only 33% of the OTU that had a relevant effect on ADG AL had a positive effect on this trait.
Additional file 1: Table S1. Posterior means, posterior medians, 95% highest posterior density intervals (HPD 95% ) structural coefficients different from zero (in g/d, CSS-normalized OTU units) that had an effect on the average daily gain of growing rabbits fed ad libitum and their assignment at the lowest taxonomic level. Table S2. Posterior means, posterior medians, 95% highest posterior density intervals (HPD 95% ) structural coefficients different from zero (in g/d, CSS-normalized OTU units) that had an effect on the average daily gain of growing rabbits under restricted feeding and their assignment at the lowest taxonomic level.